Electronic resonance states in metallic nanowires during the breaking process 
simulated with the ultimate jellium model 



E. Ogando/'l T. Torsti,2.3 N. Zabala,i'4 and M.J. Puska^ 

'Elektnka eta Elektronika Saila, UPV-EHU 644 P-K., 48080 Bilbo, Spam 
'^Laboratory of Physics, Helsinki University of Technology, P.O.Box 1100, FIN-02015 HUT, Finland 
''CSC - Scientific Computing Ltd., P.O. Box 405, FlN-02101 Espoo, Finland 
Donostia International Physics Center (DIPC) 1072 P.K., 20080 Donostia, Spain 

(Dated: February 1, 2008) 

We investigate the elongation and breaking process of metallic nanowires using the ultimate jellium 
model in self-consistent density-functional calculations of the electron structure. In this model the 
positive background charge deforms to follow the electron density and the energy minimization 
determines the shape of the system. However, we restrict the shape of the wires by assuming 
rotational invariance about the wire axis. First we study the stability of infinite wires and show 
that the quantum mechanical shell-structure stabilizes the uniform cylindrical geometry at given 
magic radii. Next, we focus on finite nanowires supported by leads modeled by freezing the shape of 
a uniform wire outside the constriction volume. We calculate the conductance during the elongation 
process using the adiabatic approximation and the WKB transmission formula. We also observe 
the correlated oscillations of the elongation force. In different stages of the elongation process two 
kinds of electronic structures appear: one with extended states throughout the wire and one with 
an atom-cluster like unit in the constriction and with well localized states. We discuss the origin of 
these structures. 



I. INTRODUCTION 

The miniaturization of the electronics components is of 
great importance in the development and improvement of 
new devices for applications in a wide number of fields. 
Although the laws of Nature are the same for macro- 
scopic and mesoscopic systems, the miniaturisation pro- 
cess is achieving the limit where the quantum behavior 
of matter starts to play an important role. 

If the size of the system under consideration is only a 
few nanometers, the atomic character of matter emerges 
and it can not be considered as a continuum. The regime 
of quantum behavior is reached also if one of the spatial 
dimensions of the system is reduced down to the Fermi 
wavelength of the conducting electrons. Then, the con- 
finement splits the continuous electronic band in this di- 
rection into a set of discrete energy levels. In both cases, 
the behavior of the system changes from what is expected 
from the macroscopic case. In metallic nanowires the 
Fermi wavelength is of the same order of magnitude as 
the atomic distance, and both atomic and electronic dis- 
crete character compete and/or couple, determining the 
properties of nanowires. 

There are many experimental and theoretical works 
which have gone deep into the understanding of the main 
features of nanowires. Experimental studies have fo- 
cused on the investigation of the mechanical and elec- 
tronic properties such as force, atomic structures and 
conductance, pointing out the close relation between 
them. Among the experimental setups we want to £. 
phasize the role of the scanning tunneling microscopeEli 
(STM) and the mechMHcally controllable break-junction 
(MCBJ) techniqucfU'EHa. In both techniques metallic 
nanowires are produced by putting two protrusions in 
contact and then pulling them away from each other over 



atomic distances. In this process, a nanowire is produced 
which upon pulling is elongated and narrowed until it 
eventually breaks. These methods have allowed the study 
of transport properties and stability of nanowires. 

The MCBJ techniques have demonstrated|-|tlie exis- 
tence of electronic and atomic shell strUjCturesclElEl, anal- 
ogous to those found in atomic clustersulj. In these ex- 
periments the conductance has been studied by building 
histograms of the conductance during the breaking pro- 
cess. The results show that there are conductance values 
which are much more probable than others. Due to the 
relation between the conductance and the radius at the 
narrowest part of the nanowire, that means that there 
are magic radii with enhanced stability while other radii 
are less stable, and therefore they appear less frequently 
in the conductance histograms. The atomic structures 
of nanowires in the last steps befo£&.|hreajdpg have been 
studied also with these techniquesElcfE30ll3. 

The experiments discussed above have been accom- 
panied by supporting theoretical investigations, which 
can be split in two groups. The first group includes 
classical and ab initio molecular dynamics simulations, 
in which the atomic structure of nanowires is taken 
into account. These investigations have been success- 
ful in many aspects, e.g., showing the atomistic mecha- 
nisms of the narrowing process (appearance of disloca- 
tions, order-disorder stages, etc.) and their link with 
other measurable quantities such as the elongation force 
or the conductanceE3EJ. Moreover, from the viewpoint 
of the present work we notice the predictions of spe- 
cial atomic .. aijrapgements in STM tips and nanowire 
neckalj'Lj'E^Oliflllj. The second group of models is more 
related to properties due to the confinement of electrons 
in reduced dimensions, and ignores the atomistic struc- 
ture of matter. In these calculations analytic approxima- 
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tions as well as self-consistent electronic-structure mod- 
els have been used, mainly within the jellium framework. 
The results obtained with these methods are also enlight- 
ening, explaining the cohesive and electronic transport 
properties of nanowires, especially in the fjf^fif [pf, alkali 
metals with strong free-electron characterElEil'LjO'Efl. 

The aim of this paper is to simulate the breaking of 
nanowires. For this purpose we choose the jellium model 
and the self-consistent electronic-structure calculations 
within the density-functional theory. In spite of their 
simplicity jellium models have provided a simple and 
transparent way to understand the physics of metallic 
nanowires. More specifically, we use the ultimate jel- 
lium (UJ)— model. This model was first proposed by 
Manninencj to investigate the structures of alkali metal 
clusters. It kaaJDcen used for the same purpose also in 
later studiesOLll. To our knowledge the present work 
is the first time the UJ-model is used to simulate the 
nanowire breaking. In practice, we solve the ensuing 
Kohn-Sham equations in a real-space pnirit grid using 
the powerful Rayleigh Quotient MultigridEaEf (RQMGl 
method implemented in the program package MIKAe2I 
(Multigrid Instead of K-spAce) . 

Within the UJ approach, the background positive 
charge density is fully relaxed in shape and density so 
that it equals at every point with the electron density. 
One can think that this freedom of the positive back- 
ground charge mimics the efficient rearrangement and 
diffusion of ions at temperatures close to the melting 
point at which the shell- and supershell-structure stud- 
ies by the MCBJ techniques have been performed for 
alkali metal^O. In principle, there is no restriction 
for the geometry of the constriction. This is in con- 
trast with the previous jellium calculations, which intro- 
duced ad hoc shapes for the nanowire. In our model the 
electrons themselves acquire self-consistently the shape, 
which minimizes the Kohn-Sham energy functional, and 
carry along the positive background charge. However, in 
order to reduce computational demands and to highlight 
the important phenomena from the complexity of possi- 
ble solutions, we restrict the shapes of nanowires to the 
axial symmetry, i.e. rotational invariance with respect 
to an axis. 

One of our main results is that in the narrowest Ptrt 
of the nanowire, electronic cluster derived structurest3't2l 
(CDS's) appear. This tendency of electrons to form em- 
bedded clusters in the jellium constrictions is analogous 
to the preferred cluster-like arrangements of atoms in 
contacts described by first-princinies atomistic calcula- 
tions by Barnett and Landmar£3M_. CDS's have later 
been reported also by other authoralj. The main differ- 
ence is that in our jellium model the atomistic charac- 
ter of the previous works is lost and the electrons alone 
are responsible for the phenomenon. The single-electron 
states provided by the jellium model can be studied in 
order to gain insight into the localization effects associ- 
ated with the CDS. The conductance of the constriction 
can be estimated either by counting the bands crossing 
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the Fermi level, or by using the WKB formula. 

The rest of the paper is organized as follows: in Section 
we describe the practical features of the UJ-model and 
the RQMG-method to calculate the electronic structure 
during the elongation process. In Section III we discuss 
the results for the electronic properties. As a starting 
point, we consider the results for infinite wires. Then we 
focus our attention on the breaking process of a finite 
cylindrical UJ nanowire supported by leads. Section IV 
contains the conclusions. 



II. THEORY 

A. Jellium models 

The jellium model has been widely used in self- 
consistent electronic-structure calculations of nanostruc- 
tures. It simplifies the problem by replacing the ions by a 
uniform rigid positive charge density background, which 
globally neutralizes the electron negative charge. The ef- 
fective potential of the Kohn-ShamEil equations is written 
as (Atomic Hartree units are used throughout this paper 
to write the equations) 

Kff = / , 71 dr + Vs,c[n-{r)\, (1) 



where the first term on the right-hand side includes 
the electron-electron and electron-positive background 
Coulomb interactions and the second term gives the 
exchange-correlation potential within the local density 
approximationESO (LDA). 

Different types of jellium approaches have been intro- 
duced. The simple jellium (SJ) model has the prob- 
lem that there is only one equilibrium charge density, 
at Ts ~ 4.18 uq (n_ = 3/(47rr^)) corresponding approxi- 
mately to the average conduction electron density in Na 
or K metals. This means that for rg values lower (higher) 
than ~4.18 ao the jellium system tends to expand (com- 
press). In the SJ-model the electron density has the same 
mean value as the positive background due to the elec- 
trostatic forces. The SJ-model gives incorrect values for 
properties such as the cohesive energy, surface energy and 
bulk modulus, due to the trend of the system to compress 
or expand. To improve, the results corrections can be 
added to the SJ-modelH, e.g. using the jSo-called stabi- 
lized jellium model introduced by PerdewE3 and ShoretJ. 

In this work we use the UJ-model, the philosophy of 
which differs from the stabilized jellium model in that 
it does not try to correct the above-mentioned deficien- 
cies of the SJ-model. The peculiarity of the UJ-model is 
that the positive charge background is allowed to relax. 
The UJ-model represents the ultimate limit in which the 
positive background is completely deformed to have the 
same density as the electrons locally at every point. In 
this way, the Coulomb term in the potential always van- 
ishes, and in Eq. (1) only the exchange-correlation term 
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survives. The total energy is then minimized in the in- 
terplay between the exchange-correlation and the kinetic 
energies. 

One limitation of the UJ-modcl is that, as in the SJ- 
modcl, there is only one equilibrium charge density, at 
rs ~ 4.18 flQ. But, the absence of electrostatic potential 
disables the mechanism to keep the electrons at a given 
density, and inside the UJ the mean electron density be- 
comes equal to the equilibrium density. Another property 
of the UJ-model, derived also from the absence of electro- 
static potential, is that the shape of the electron density 
is to a large extent uncontrollable, and it evolves until the 
ground state is achieved. This property has been used to 
study the. mo.stj favorable shapes of simple- metal atom 
clustersE^E^LllcZl. In the present work, however, we have 
to deal with open systems and we have to impose certain 
controlling restrictions in order to model the pulling of 
the nanowires. The description of the solutions to these 
requirements is postponed to Section [II C. 



along the wire axis, the relaxation of the positive back- 
ground charge and electron density is limited in the radial 
direction. Consequently, it is necessary to solve numer- 
ically only the rad-ial part of the Schrodinger equation 
(see Zabala et alx3 for technical details). 

For the systems studied in Sections III B| and III C , 
however, the translational invariance is not required. 
But, in addition to the rotational invariance, periodic- 
ity in the axial direction is assumed with unit cell length 
Lcoii. Thus, the wave functions ^p are indexed by the 
quantum numbers m, n, and kz- Here, m is the angular 
momentum quantum number and kz is the Bloch wave 
vector along the wire axis. With m and kz given, n enu- 
merates the orthogonal states in the order of increasing 
energy eigenvalue. The UJ system is solved by finding 
the self-consistent solution to the set of equations 
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B. Numerical methods 



In Section HI A infinite uniform cylindrical wires are 



studied. Since these systems are translationally invariant 



Unik^nir, z + Lcoii) = e*'=--^""[/„fc^„(r, z). 
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n(r) = 2 ^ (2 - 5om)/mfe.„|C/mfe.„(r)P, (5) 

mkzTL 

V^?i(r, z) = Kclr, z) = — — ^ — . (6) 

on(r, z) 

The effective potential Vcff (r,z) equals the exchange- 
correlation potential V^c(^, The electron density 
n(r, z) is obtained by summing single-electron densities 
with the occupation numbers fmk^n- The degeneracies of 
the states are taken into account by the factor 2(2 — Som) 
and the occupation numbers fmk^n obey the Fermi-Dirac 
statistics with the Fermi level {Ep) so that the system 
is neutral. A finite temperature of 1200 K is used to 
stabilize the solution of the set of equations. 

The Schrodinger equation (^ is discretized on a reg- 
ular two-dimensional (r, z) point mesh. We use stan- 
dard fourth-order central-difference discretizations for 
the first and the second derivatives. The grid is sur- 
rounded by a frame with the thickness of two grid points. 
These ghost points are necessary for the evaluation of 
the derivatives near the edges of the computation vol- 
ume. The wave functions are required to vanish at the 
ghost points corresponding to the radial surface of the 



cylindrical computation volume, whereas at the axis the 
values at ghost points can be evaluated by noting that 
U{-r, z) = (-l)'"J7(r, z). The Bloch-condition (Eq. (|)) 
gives the recipe for obtaining the values at the ghost 
points of the periodic boundary. 

The problem with standard real-space relaxation 
methods for Eq. (^ is the so called critical slowing- 
down phenomenon resulting from the fact that at a time 
they use information from a rather localized region of 
space. As a result of the locality, the high-frequency er- 
ror, corresponding to the length scale of the grid spacing, 
is reduced very rapidly in the relaxation. However, once 
the high frequency error has been effectively removed, 
the very slow convergence of the low-frequency compo- 
nents dominates the overall error reduction rate, i.e. crit- 
ical slowing-down occurs. Multigrid methods avoid this 
problem by treating the low frequency components of the 
error on coarser grids, where their wavelength is compa- 
rable to the grid spacing. 

Applying the multigrid methods to the Schrodinger 
equation is a fairly complicated task because one has to 
solve both the eigenvalue and the wave function simul- 
taneously - this makes the problem nonlinear. Also, one 
has to solve several wave functions simultaneously, avoid- 
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ing the bottleneck of orthogonalizations as well as possi- 
ble. Staadard methods based on the fuU-approximation- 
storageES method require, that the wave functions are 
well representable on the coarsest grid used, implying 
severe limitations on the acceleration obtained by the 
multigrid idea. We use the recently developed gen- 
eralizatiojx-pf the Rayleigh-quotient multigrid (RQMGj) 
methodtElo as implemented in the MIKA-packageLJ, 
that avoids the problems described above. In short, one 
applies the Gauss-Seidel method on the finest grid. On 
the coarser grids one applies coordinate relaxations on 
the functional 



(V'nlV'n) 



i(V'.i^n)r 



(7) 



This functional, which is actually defined on the finest 
grid, is the sum of the Rayleigh quotient and a penalty 
functional, which is introduced to insure the orthogonal- 
ity. Moreover, the relaxations are performed sirnultane- 
ously for all wave functions. See Heiskanen et alsn for a 
more thorough discussion of technical details. 

The Kohn-Sham equations have to be solved self- 
consistently. In other words one has to iterate until the 
output potential Vos obtained from Eq. (|^) equals the in- 
put potential Ves that is used in Eq. (||). In typical cases 
of electronic structure calculations, to avoid divergence 
due to charge sloshing, one uses sophisticated strategies 
to construct the input potential for the next iteration as 
an optimized mixtur.e of input and output potentials of 
previous iterationsOCj. In the UJ-iterations, however, 
the output potential can be taken directly as the input 
potential of the next iteration resulting in a rapid conver- 
gence. This is because of the absense of the long-range 
Coulomb interaction, which is the cause of the charge 
sloshing phenomenon. 



III. RESULTS AND DISCUSSION 

A. Infinite uniform cylindrical wires 

The main results of this paper concerni ng the nanowire 
breaking process are discussed in Sections III B and III C . 
As a preliminary work, and in order to gain insight into 
the UJ-modcl in comparison with the stabilized jellium 
model, we study the stability of infinite uniform cylindri- 
cal UJ- wires. 

We calculate the surface energy of the nanowires and 
the oscillations in the energy per unit length as it was 
made in our previous work describing Al, Na-^apd Cs 
nanowires within the stabilized jellium modcltjO. The 
results are shown in Fig. |l| as a function of the wire 
radius R. Here the radius is defined as the radius of 
the positive background charge in the SJ system with 
Ts = 4.18 ao and the same amount of charge per unit 
length. In order to separate the energy oscillation from 
the average behavior, the so-called liquid-drop modclcll 
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FIG. 1: Surface energy of infinite uniform cylindrical UJ 
nanowires as a function of the nominal wire radius (See text). 
In the inset the energy oscillations are shown and the first 
magic radii are marked. 



is used. In this model the energy of the jellium system 
can be written as the sum of two terms - one propor- 
tional to the volume and the other proportional to the 
surface area. For the first term, the energy/volume ra- 
tio -corresponds naturally to the homogeneous electron 
gascj with — 4.18 (ta-j-Jiis view has been tested in 
clusters and nanowire^3E2li23 and it describes correctly 
the mean energy, i.e. without the characteristic oscilla- 
tions due to the quantum confinement. We fit the self- 
consistently calculated total energy per unit length to a 
liquid-drop model type function. Then, subtracting this 
smooth energy function from the total energy we get the 
pure energy oscillations, which are shown in the inset of 
Fig. 0. Note that there are radii for which the energy is 
at minimum. They correspond to wires which are more 
stable than wires with slightly different radii and higher 
energies. The first magic radii are at i? = 4.3, 7.3, 10.3, 
13.6, 17.8, and 20.7 oq. We use these radii for the ini- 
tial uniform wires in the nanowire breaking simulation in 
Section IIIC. The shell and-super-shell structures stud- 
ied in previous calculationsEJO are also quite clear. In 
comparison with the energy oscillations of Na we observe 
that the beat positions are shifted to higher radii. The 
reason is that the UJ potential is softer at, the surface 
than the stabilized jellium potential for Nan3. 



B. Periodic systems 

Now we change the scheme and allow the wire to de- 
form also in the axial direction. However, we impose pe- 
riodic boundary conditions with the unit cell length LccW 
along the wire axis. Lcdi is thus the maximum pertur- 
bation wavelength in our calculation. From the liquid- 
drop model point of view, neglecting the small contri- 
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FIG. 2: Periodic infinite UJ wire with the nominal radius of 
R = 5.5 ao and 6 (a), 7 (b), 8 (c), and 10 (d) electrons in the 
periodic unit cell. The figures show the electron density of 
one unit cell. The contour spacing is 0.15 times the UJ bulk 
density value. 



bution of the curvature energy, the liquid wire attempts 
to achieve the shape which minimizes the surface, and 
thereby the total energy. Under this assumption an infi- 
nite periodic liquid wire is a uniform cylinder for lengths 
LccW < 4.5i?. For Lceii > 4.5ii! it deforms trying to 
achieve the energetically most favorable state, an infinite 
chain of spheres. However, Kassubek et alE3 showed us- 
ing a semi-classical model and perturbation theory that 
due to the discreteness of electronic structure the wires 
with magic radii remain uniform also at large Lccu/ R ra- 
tios. With this result they argued that in the narrowing 
process of an infinite wire, when the radius is crossing an 
unstable zone before the next stable radius is achieved, 
the wire would spontaneously deform acquiring a wavy 
or deformed shape. We corroborate these results non- 
perturbatively using the UJ model as follows. 

We choose a certain radius R and solve for the UJ elec- 
tronic (and positive charge density) structures imposing 
increasingly longer supercell lengths Lceii by increasing 
the number of electrons in the cell. Thereby we deter- 
mine the critical supercell length (the wavelength of a 
perturbation) at which the wire starts deforming. For 
magic wires we find no wavy solutions, the wires remain 
uniform. For example, for R = 7.3 oq the wire is still 
uniform at Lcdi/R ~ 36. The wires corresponding to 
the radii at the maxima of the energy oscillations in Fig. 
|l| are the most unstable ones. These wires are uniform 
up to a critical value of Lcdi but above it they sponta- 
neously deform to a wavy or non-uniform density profile 
along the wire axis. As an example. Fig. |2| shows the 
behaviour of a wire with radius R — 5.5 ao when the 
number of electrons in the unit cell is 6, 7, 8, and 10 
and the unit cell lenght Lcoii increases as 19.3, 22.5, 25.7 
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FIG. 3: Schematic view of the model system for simulations 
of breaking of finite nanowires supported by two leads. 



and 32.1 ao, respectively. The unit cell with eight elec- 
trons correspond to a magic spherical cluster and that 
of ten electrons the pair of magic clusters of eight and 
two electrons. The critical values for the unstable radii 
of i? = 5.5, 8.6, 11.6 and 19 ao are Lcdi/R =4.1, 3.2, 4.2, 
and 4.8, respectively. I.e., we obtain values near the clas- 
sical value of 4.5. At the unstable radius of i? = 15.5 ao 
the wire is not deformed at least up to Lcdi/R = 10 (the 
largest length we have calculated), probably due to the 
fact that this radius lies in a beat of the super-shell struc- 
ture and it is actually relatively stable. We start all the 
calculations with a converged uniform potential profile 
along the wire axis (See Section III A ). In this way we 
do not "add any energy" to the system when initiating 
the calculation. Therefore, if the wire starts to deform in 
the iteration process the reason is the disappearance of 
the local energy minimum. 

In addition, we narrow a stable uniform wire by in- 
creasing the length Lcdi of the periodic cell and maintain 
the number of electrons constant. Each elongation step 
is solved self-consistently until convergence is reached. 
We observe that during the first steps the wire remains 
uniform but at some point, before breaking into isolated 
clusters, the wire spontaneously deforms. Thus, we con- 
firm self-consistently and dynamically the hypothesis by 
Kassubek et alE3. 



C. Breaking of supported finite nanowires 

In order to study the formation and evolution of 
nanoconstrictions between two supporting leads we fol- 
low the next procedure. First, we fix the number of elec- 
trons in the periodic supercell and solve self-consistently 
for the electronic structure of a uniform UJ wire having a 
stable magic radius. Then, the potential at both ends of 
the periodic cell is "frozen". This means that, although 
the Kohn-Sham equations are solved in the whole wire, 
in these regions the potential is not updated in the self- 
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FIG. 4: Supported UJ wire. The UJ constriction contains 
eight electrons. Density contour plots for four different elon- 
gation lengths: AL — 7.9 ao (a), 19.8 ao (b), 20.8 ao (c), and 
25.8 ao (d) are shown. The snapshots in (b) and (c) are from 
consecutive self-consistent calculations and the snapshot (d) 
is the last step before the nanowire breaking. The contour 
spacing is 0.15 times the mean UJ bulk density value. 



consistency process. The function of this "frozen" part is 
to emulate the lead parts where ion rearrangement does 
not occur as efficiently as at the constriction. In our cal- 
culation, these leads serve as handles to grab the UJ and 
pull it. The rest of the wire, the UJ at the middle part 
of the supercell, is the place where the wire will stretch. 
A sketch of the configuration is shown in Fig. ^. A sharp 
change in the potential between the constriction and the 
leads turned out to cause difficulties in numerical calcu- 
lations. Therefore, we smooth out the potential at the 
left edge using the form 

~ 4dgo)^frozcn + ^^C^dgo - z)Vuj, (8) 

where F is a Fermi function with half- with of 0.5 ao, 
^frozen is the " frozcu" potential and Vuj is the self- 
consistent UJ potential. For the right edge an analogous 
mixing is used. The main properties of the nanowire 
will not depend on the particular choice of this match- 
ing because the physical features are determined by the 
narrowest part of the constriction. 

We perform simulations starting with radii between 
7.3 flo and 20.7 ag and changing the number of electrons 
initially in the constriction. The elongation of the wire 
is made in steps of about one atomic unit, and always 
starting from the previous converged density, so that the 
grid spacing of the point mesh is increased to enlarge 
the cell. In order to overcome the intemctions between 
the constriction and its periodic replicacj, we choose the 
length of the lead part to be 6 Fermi wavelengths {Xf = 
13.7 ao). Throughout the rest of the paper we will use 
AL for the elongation; AL = for the first step. 



In Fig. we show snapshots of the electronic density 
for a wire with the starting radius of 10.7 Oq. The UJ- 
part corresponds to eight UJ-electrons placed initially in 
the neck region. Electrons are free to move inside or 
outside the leads, depending on the requirements of the 
self-consistent solution. However, there are always about 
eight electrons in the constriction. Although this is one 
of the smallest wires we have calculated, it shows all the 
main features observed when simulating also larger wires. 

If the breaking of an UJ nanowire would happen as 
for fluid between the leads, the electron density should 
evolve forming sucatenoid-shaped surface. Similar shapes 
(like hyperbolicE^, paraboliccj, cosinal3, etc.) have been 
used before to model the nanoconstriction in simple free- 
electron or jellium simulations. The main results, when 
the comparison is possible, have been essentially the same 
irrespective of the actual shape. In Fig. ^(a) the electron 
density is shown after the elongation of AL — 7.9 ao. 
The catenoid-like density profile appears as expected 
for a classical fluid. When we continue elongating the 
nanowire the shape of the electron density changes dra- 
matically from the classical one. If the distance between 
the leads is short the electrons are strongly trapped at 
the narrowest part and they do not have much freedom in 
the rearrangement process. When the length of the con- 
striction is large enough, the electrons have more space 
and freedom to achieve different types of energetically 
preferred shapes. 

In Fig. ||(b) AL = 19.8 Oq and the electrons in the 
constriction form a CDS. The electron density per unit 
length has two minima at both sides of the CDS and there 
are 7.1 electrons between these narrowest cross sections. 
The embedded cluster reminds the closed-shell cluster of 
eight electrons, but there are some differences. There are 
not enough electrons and the symmetry is not exactly 
spherical. It seems that the pz orbital (z along the cylin- 
der axis) of the cluster has disappeared. We will analyse 
the structure in more detail below. Figure ^(c) shows 
the next consecutive elongation step with AL = 20.8 ao. 
Note that the CDS disappears and a sudden change in the 
mean radius happens. In fact, the conductance changes 
simultaneously abruptly from 3 Go to 1 Go (See the in- 
set in Fig. ||(a)). At this point it is also remarkable 
that the shape of the constriction is again far from the 
catenoid having a constant magic radius. Figure ^(d) 
is for AL = 25.8 ao, the last step before the nanowire 
breaks. Again a CDS appears during the elongation from 
the third to the fourth snapshot. There are 1.8 electrons 
between the two minimum cross sections at both sides of 
the CDS. This CDS can be interpreted as an embedded 
two-electron cluster. We observe that the radius of the 
constrictions is more or less constant with the same value 
as in the previous snapshot in Fig. ^(c). 

At this point we want to focus on one characteristic 
property of UJ found when simulating the wire breaking: 
the UJ-matter deforms very easily. This ability to de- 
form allows the formation of the cylinders of magic radii 
glued to the leads. The radius jumps from one magic 
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FIG. 5: Main figures: conductance, effective radius at the 
constriction, and elongation force for a wire with initial radius 
R = 20.7 ao and about 60 UJ-electrons in the constriction. 
Insets: the same quantities for the wire in Fig. ^ with initial 
radius R — 10 ao and 8 UJ-electrons in the constriction. The 
arrows mark the points were the density has been plotted in 
Fig. I 



radius to the next through an abrupt charge reorganiza- 
tion. The CDS's of about two or eight electrons appear 
before the last charge reorganizations and the wire break- 
ing. If there is enough UJ between the leads suspended 
long thin cylinders appear and in the last steps they al- 
ternate with chains of CDS's producing a very extended 
elongation process. Here we want to underline that the 
CDS formation is a process different from the stability 
of a uniform cylindrical wire against the formation of a 
chain of spheres studied in the previous Section IIIB. 



In that section, the quantum-mechanical shell structure 
may conserve the cylindrical structure which is not classi- 
cally stable whereas now the quantum-mechanical shell- 
structure effect destroys the classical catenoid-type of so- 
lution producing a CDS in the constriction. 

In Fig. ||we show the conductance, the effective radius, 
and the elongation force as a function of the elongation 
for two different wires. The main figures correspond to 
an initial configuration with the radius of 20.7 gq and 60 
electrons in the UJ-constriction. The insets display the 
results for a wire with an initial radius of 10.7 oq and 
eight electrons in the constriction. The electron density 
of the latter wire is plotted in Fig. ^ at certain elongation 
stages. 

The conductance is calculated with the adiabatic a: 
semiclassical approximation used by Brandbyge et al. 
The constriction is divided into transversal slices. Then 
for each slice a uniform wire with the radial extent of the 
slice is built and the energy eigenvalues of the subband 
bottoms are calculated for this slice. The subband bot- 
toms give effective potentials along the wire axis. If we 
look at the dependence of one of them on the position 
we see that it raises at the constriction due to the strong 
confinement (see Fig. The electrons in this subband 
at the Fermi energy of the leads have to overcome this 
barrier in order to carry current. To evaluate the trans- 
mission probability of the electrons at the Fermi level 
through the barrier the semi-classical WKB formula is 
used. 

The properties of the nanowires have been demon- 
strated to be dominated by the narrowest part of the 
constriction. Therefore we calculate an effective radius 
by evaluating the electron density per unit length at the 
middle of the wire. It is obtained with the value of the 
bulk electron density (corresponding to = 4.176 ao). 
Fig. ||(b) shows the effective radius as a function of the 
elongation of the wire. The plateaus or shoulders are 
in good coincidence with the infinite wire magic radii of 
10.3 flo, 7.3 ao and 4.3 ao. For the larger wire shown in 
Fig. |(b) also a small kink can be seen at AL = 15.5 ao, 
which corresponds to the magic radius of 13.6 ao. Wider 
magic radii do not appear because of the beat region of 
the supercell structure. In the inset at the end of the 
plateaus the effective radius increases when elongating 
the wire due to the CDS formation. The sudden decrease 
of the effective radius, accompanied by a step in the con- 
ductance, is due to the sharp charge rearrangements in 
the constriction. 

The elongation force, shown in Fig ||(c) is evaluated 
as the negative derivative of the total energy with re- 
spect to the elongation. The rearrangement of the wire 
charge leads to discontinuous upward steps in the force, 
while if the radius changes smoothly the force draws a 
continuous buckling curve. At this point we want to 
point out the superiority of the UJ jmodeL in the force 
calculation over 0|t|ljier jellium models£30'0. In contrast 
with experimentsErH, the latter show a continuous be- 
havior of the force without any steps. Moreover, for 
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FIG. 6: (a) Electron band structure of the wire having eight UJ-electrons and the elongation of AL = 19.8 ao (Fig. ^(b)). The 
vertical dashed lines mark the different Brillouin zones. The label for each branch represent the (m, n) subband for the infinite 
wire. The energy eigenvalues are solved at two k2-points: at the origin and at the zone boundary of the supercell Brillouin 
zone, (b) LDOS integrated between the two narrowest points of the electron density in Fig. ^(b) is displayed. The solid line 
represents the total ILDOS, the dashed lines are the contributions to the ILDOS of the states with m=0 (long-dashed line) 
and m=l (short-dashed line). The dotted lines show the DOS of hypothetical localized states. The states marked with squares 
in the band structure are the main states contributing to the ILDOS. In the inset, the analogous plot for the ILDOS in the 
constriction of Fig. ^(a) with AL = 7.9 ao is shown. The origin of energy is the Fermi level. 



narrow constrictions positive values are obtained when 
the wire crosses an unstable zone. Note that in our 
model the fefce is always negative, as obaeswedjn the 
experimentaJ'cl and in atomistic simulationgjllj'L3'LZl. Fig. 
H shows clearly that the transport, geometrical and me- 
chanical properties of the nanowires under elongation are 
related. 



D. Electronic cluster-derived structures 

Let us now analyse more closely the CDS appearing 
in Fig. ^(b). In order to enlighten the origin of this 
structure we plot in Fig. ^a) the single-particle energy 
spectrum of the wire. The extended zone scheme is used 
for clarity. The labels on the left of each branch repre- 
sent the corresponding (|m|,7i) subbands for the infinite 
wire. In practice, n is obtained by calculating the num- 
ber of radial nodes at the cell boundaries (See Fig. ||). 
The branches have the characteristic parabolic shapes, 
but they show two different stages. In the lower part of 
the parabolic subbands the eigenvalues form flat plateaus 
without kz dispersion. These states correspond to wave- 
functions localized at the leads and they vanish at the 
center of the constriction. Therefore the (0,2) and (2,1) 
subbands can not carry current through the constriction 
and they are closed channels. On the other hand, the 
states of the upper part of the (0,1) and (1,1) branches 
are extended along the whole wire and they form a con- 
tinuous band (with the exception of small band gaps). 



The conductance of the wire is thus 3 Go due to the ex- 
tended states of the (0,1) and (1,1) open channels at 
Fermi energy. This conclusion is in accordance with the 
value obtained with the WKB approach. 

In Fig. ^(b) we plot the integrated local density of 
states (ILDOS) in the constriction for the band struc- 
ture of Fig. ^(a). It is calculated by integrating the local 
density of states (LDOS) over the space between the two 
narrowest parts in the electron density in Fig. ^b) . The 
LDOS itself is obtained by substituting the discrete en- 
ergy levels with Lorenzians of width (FWHM) of 0.4 eV 
and weighting them by the local probability amplitudes 
of the states in question. The ILDOS has two clear peaks, 
and when decomposing it we can see that the lower and 
the higher peak have the m = and m = 1 character, 
respectively. The contribution of the m = 2 states is 
negligible. The two ILDOS peaks can be fitted by two 
energy levels convoluted with the same Lorenzian as the 
eigenlevels in the LDOS calculation. The resulting reso- 
nance peaks are shown in Fig. ^(b) by dotted lines. The 
positions and the heights of these peaks have been fitted 
manually. The coincidence between the fit and the true 
ILDOS is remarkable. In the inset of Fig. ^(b) we plot 
the LDOS integrated between the leads for the electron 
density showed in Fig. ^(a) having no CDS. We observe 
that the ILDOS is much smoother and it is similar to 
the DOS of an infinite wire with delocalized states. The 
different m contributions cannot be fitted by single reso- 
nance peaks as shown by the dotted peak for the m — 
and m = 1 contributions. Moreover, the inset shows that 
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FIG. 7: (a), (b) and (c): Selected single-electron states in the wire having eight UJ-electrons and the elongation of AL — 19.8 ao 
(Figs, ^(b) and Contour and profile plots (along the wire through the maximum value) of the density are shown. Plot 
(a) corresponds to a m = resonance state at the low energy peak in the ILDOS (Fig. Plot (b) corresponds to a m = 1 
resonance state at high energy. Plot (c) is an extended state of the m = subband at the energy of the m = 1 peak in the 
ILDOS. Plot (d) is a localized state with m = corresponding to the elongation of AL = 25.8 ao (Fig. §(c)). The contour 
spacing is one tenth of the maximum value. The kz vector is given in Brillouin zone units (-^). 



the m-decomposed peaks are slightly asymmetric with a 
tail on the high energy side. These tails, which are not 
observable in the main figure in which the CDS appears, 
are due to the y^-dependence of the subband peaks in 
the DOS for infinite wires. 

The ILDOS analysis suggests that in the energy sub- 
bands or branches, at the transition points from states 
localized in the leads to states extended across the whole 
wire (see Fig. §(a)), rather localized resonance states ap- 
pear in the constriction. To clarify this point we plot 
selected states at the ILDOS peak energies in Fig. 0. 
Figs. |^(a) and |^(b) show clearly the localized character 
of the wavefunctions in the constriction at these energies. 
The state in Fig. ^(a) can be identified as the Is orbital 
of an eight-electron cluster. The second well-localized 
state, Fig. ^(b), has m = 1. Therefore it is doubly de- 



generate and it is identified as the p^y orbital. At about 
the energy of this p^y orbital a p^ orbital (directed along 
the wire axis) should appear in the m = branch in 
order to complete the eight-electron cluster. However, 
we do not find such a state with a strong localization in 
the constriction. As shown in Fig. ^(c) the pz type states 
are much more delocalized than the pxy resonance states. 
The difference reflects the fact that due to the orienta- 
tion the interaction of the cluster pz orbital with the lead 
states is much stronger than that of the p^y orbital. The 
absence of a well-localized Pz orbital explains the clearly 
non-spherical shape of the embedded cluster in the elec- 
tron density plot of Fig. |^(b), and also the finding that 
there are only 7.1 electrons in the constriction between 
the two narrowest cross sections. 

Fig. |^(d) shows a well-localized state for the wire with 
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FIG. 8: Effective ^-dependent potentials for the different 
(m, n) channels. The wire has eight UJ-electrons and the 
elongation of AL = 19.8 ag (Figs, ^(b) and |). The dashed 
lines in the potential wells are drawn at the energies of the 
resonance states. The lowest line is the bare Kohn-Sham po- 
tential at = 0. 



eight UJ-electrons at the elongation of AL = 25.8 oq 
(Fig 1(d)). The state can be identified as the Is orbital 
of a 2 electron cluster glued to the leads. There are 1.8 
electrons between the two narrowest cross sections of the 
constriction supporting the assumption that this state is 
related to a two-electron cluster. 

The states in Fig. ^ have always a wavy, non-decaying 
background. This is a characteristic of resonance states; 
truly localized states would decay exponentially. The 
wavy background corresponds to the wave function of 
the leads (plane wave) at the energy, which matches with 
that of the cluster state. To check this assumption we re- 
alize that the wavelength of the plane-wave background 
corresponds to the kz quantum number in the extended 
zone scheme; There are indeed two maxima in the mod- 
ulus of the wave function per every Brillouin Zone unit 
of kz (see the labels of each wave function) . 

The existence of resonance states is related to the 
shape of the self-consistent potential having a small po- 
tential well in the nanoconstriction. To point out how 
this potential can admit a resonance state we show in Fig. 
^ the effective potential for states with different (m, n) 
quantum numbers, calculated within the adiabatic ap- 
proximation for the wire with eight UJ-electrons and the 
elongation of AL = 19.8 aq (Fig. |(b)). We see that 
electrons at the constriction feel the existence of a po- 
tential well. We plot the energies corresponding to the 
ILDOS peaks with dashed lines and note that they lie 
exactly in the potential wells, where the resonances situ- 
ate. It is also evident that an occupied resonant pz state 
does not occur because its energy eigenvalue should be 
well above the effective potential of the (0, 1) branch and 
because the potential well of the (0, 2) branch is above 
the Fermi level. In addition, by the help of Fig. ^ we can 
explain the different parts of the electron energy bands 



in Fig. ^(a). The states with energies above the effec- 
tive potential maxima are extended along the whole wire. 
These are the current-carrying states of each branch. The 
states below the potential minimum of the constriction 
are trapped in the leads, corresponding to flat plateaus in 
the lower part of the energy branches (Fig. ^(a)). Finally, 
between the potential maxima and the local minimum in 
the center we find resonant states which are enhanced at 
the constriction although they continue as plane waves 
in the leads. 



IV. CONCLUSIONS 

We studied the stability of nanowires and the nanowire 
breaking process performing self-consistent calculations 
within the ultimate jellium model. In the model, elec- 
trons and positive background charge acquire the optimal 
density minimizing the total energy. The model enables 
thus studies of shape-dependent properties of nanoscopic 
systems such as quantum dots or, as in the present work, 
quantum wires. The model advocates the idea that the 
electronic structure determines via the shell structure the 
geometry and ionic structure also in a partially confined 
system. 

First we analysed the stability of infinite periodic quan- 
tum wires pointing out the ability of the electronic band 
structure to stabilize the nanowires at magic radii, i.e. 
any small deformation of the nanowire along the z axis 
always increases the energy. At the unstable radii cor- 
responding to maximum values of the energy oscillations 
the wire is uniform up to a critical value of the unit cell 
length. The critical values found are close to the classical 
value of Lccw/R = 4.5. Above this limit the local energy 
minimum disappears and a deformation of the wire low- 
ers the total energy. 

Then we investigated the elongation process of finite 
nanowires supported by leads. The elongation force, con- 
ductance and effective radius of the constriction were cal- 
culated simultaneously. The importance of the charge 
relaxation, in order to obtain results in agreement with 
experiments, was shown, e.g., in the case of the elonga- 
tion force. The ability of the ultimate jellium (electron 
density) to acquire the optimal shape allows the forma- 
tion of CDS's showing the importance of electron states 
in the formation of these structures. The related reso- 
nance states and their origin was also shown. We found 
CDS's that can be linked with the eight- and two-electron 
free-standing clusters. 
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